function res = trainModel_libsvm(stego, cover, ntrain, ntest, loop, verbose, c, g, ncg)
% TRAINMODEL train SVM model using LIBSVM..
% Usage:
%   RES = TRAINMODEL_LIBSVM(STEGO, COVER, NTRAIN, NTEST[, LOOP, VERBOSE, C, G, NCG])
%   This function splits training sets and testing sets, trains the model, 
% and return the testing result.
% 
% Input:
%   STEGO  : a L*N matrix for features of stego images.
%   COVER  : a M*N matrix for features of cover images.
%   NTRAIN : size of training set. (one pair of stego/cover images is 2!)
%   NTEST  : size of testing set.
%   LOOP   : number of loops. Default value is 10.
%   VERBOSE: 1 to display a lot of messages.
%   C, G   : Value for LIBSVM. Ignore it, value 'automatic' or NaN tells 
%            the program to find  best values for C,G using search
%            (and this process is very slow).
%   NCG    : size of training sample when searching for C, G. Default value
%            is NTRAIN.
% 
% Output:
%   RES : A struct with the following fields:
%       MEAN_ERR  : Mean of error.
%
res = [];
% Linear
cg_opt = '-t 0 -h 1 ';
train_opt = '-q -t 0 -h 1 ';
% RBF
%cg_opt = '-h 1 ';
%train_opt = '-q -h 1 ';

predict_opt = '';
res.CG_OPT = cg_opt;
res.TRAIN_OPT = train_opt;
res.PREDICT_OPT = predict_opt;

if nargin < 6
    verbose = 0;
end
if nargin < 9
    ncg = ntrain;
end
%% load feature matrix
if ischar(stego)
    if verbose
        fprintf('loading stego from file %s.\n', stego);
    end
    load stego -mat res
    stego = res;
    clear res;
end
[L,N] = size(stego);
if ischar(cover)
    if verbose
        fprintf('loading cover from file %s.\n', cover);
    end
    load cover -mat res
    cover = res;
    clear res;
end
[M,N1] = size(cover);
if N ~= N1
    error('trainModel:featureNumNotMatch','Number of features in STEGO and COVER does not match.');
end
clear N1
samples = [stego;cover]; % 1..L: Stego,  L+1..L+M: Cover

%% initialize output value
res.stego_sample = L;
res.cover_sample = M;
res.feature_num = N;
nSample = L+M;

% stego(1) | cover(0)
if (nargin < 8) || strcmpi(c, 'automatic') || strcmpi(g, 'automatic') || isnan(c) || isnan(g)
    fprintf('finding best c and g\n');
    if nargin < 9
        ncg = ntrain;
    end
    nt = floor(ncg/2);
    p1 = randperm(L);
    p2 = randperm(M);
    p1 = p1(1:nt);
    p2 = p2(1:nt);
    train = [stego(p1, :); cover(p2, :)];
    ltrain = [ones(nt,1);zeros(nt,1)];
    if verbose
        fprintf('train set: stego = cover = %d\n', nt);
    end
    tic;
    [acc, c, g, cvmat] = SVMcg(ltrain, train, cg_opt);
    if verbose
        disp(cvmat);
    end
    T=toc;
    fprintf('c=%g, g=%g\n', c, g);
    res.CG_ACC = acc;
    res.CG_TIME=T;
    clear train ltrain
end
clear stego cover

if (ntrain <= 0) || (ntest <= 0) || (ntrain + ntest > nSample)
    error('trainModel:invalidInput', 'NTRAIN + NTEST must not be bigger than total sample numbers.');
end
labels = [ones(L,1);zeros(M,1)];
res.NTRAIN = ntrain;
res.NTEST = ntest;

if nargin < 5
    loop = 10;
elseif loop <= 0
    error('trainModel:invalidInput','Loop must be >0.');
end
res.loop_num = loop;
errors_fa = nan(loop, 1); %false alarm 0 said to be 1
errors_md = nan(loop, 1); %miss detect 1 said to be 0
accs = nan(loop, 3);
errors = nan(loop, 1);
res.MEAN_ERR = NaN;
res.mean_fa = NaN;
res.mean_md = NaN;
res.MIN_ERR = NaN;

if verbose
    fprintf('trn/tst = %d/%d, loop = %d\n',ntrain, ntest, loop);
end
res.C = c;
res.G = g;
cloop = 0; % Count of times Training succeed

if verbose
    fprintf('train / test\n');
end

tic;
%% train and test
% it seems libsvm is not thread safe
for i = 1:loop
    perm = randperm(nSample);
    if verbose
        fprintf('Round %d/%d\n',i,loop);
    else
        fprintf('.');
    end
    % 1 .. ntrain              => train
    % ntrain+1 .. ntrain+ntest => test
    % ntrain+ntest+1 ...       => skipped
    train = samples(perm(1:ntrain), :);
    ltrain = labels(perm(1:ntrain));
    exportSVM(['train_',num2str(i),'.dat'],ltrain,train);
    if verbose
        fprintf('train set: stego/cover = %d/%d\n', sum(ltrain), ntrain - sum(ltrain));
    end
    try
        model = libsvmtrain(ltrain, train, [train_opt, ' -c ',num2str(c),' -g ',num2str(g)]);
    catch ME
        ME
        continue
    end
    clear train ltrain
    
    test = samples(perm(ntrain+1:ntrain+ntest),:);
    ltest = labels(perm(ntrain+1:ntrain+ntest));
    exportSVM(['test_',num2str(i),'.dat'],ltest,test);
    if verbose
        fprintf('test set : stego/cover = %d/%d\n', sum(ltest), ntest - sum(ltest));
    end
    [tr, acc, ~] = libsvmpredict(ltest, test, model, predict_opt);
    clear test
    
    pfa = sum((tr == 1) & (ltest == 0)) / sum(ltest == 0);
    pmd = sum((tr == 0) & (ltest == 1)) / sum(ltest == 1);
    err = (pmd + pfa)/2;
    errors_fa(i) = pfa;
    errors_md(i) = pmd;
    errors(i) = err;
    accs(i,:) = acc';
    if verbose
        fprintf('err = %0.3f\n',err);
    end
    cloop = cloop + 1;
end
fprintf('\n');
T = toc;
res.TIME = T;
res.errors_fa = errors_fa;
res.errors_md = errors_md;
res.errors = errors;
res.accuracy = accs;
if cloop > 0
    res.MIN_ERR = min(errors);
    errors(isnan(errors)) = 0;
    errors_fa(isnan(errors_fa)) = 0;
    errors_md(isnan(errors_md)) = 0;
    res.MEAN_ERR = sum(errors) / cloop;
    res.mean_fa = sum(errors_fa) / cloop;
    res.mean_md = sum(errors_md) / cloop;
end
if verbose
    disp(res)
    fprintf('Mean Accuracy: %0.2f%%\n',(1-res.MEAN_ERR)*100);
end 
end

function exportSVM(fname, lbl, inst)
    fid = fopen(fname, 'w');
    if (fid == -1)
        error('Cannot open file %s for write', fname);
    end
    [M,N]= size(inst);
    for i = 1:M
        fprintf(fid,'%d',lbl(i));
        for j = 1:N
            fprintf(fid,' %d:%0.4f',j,inst(i,j));
        end
        fprintf(fid,'\n');
    end
    fclose(fid);
end